Critical behavior of the two dimensional 2A 3A, AA binary system 



O 

o 

(N 

Ph. 
< 



Geza Odor 

Research Institute for Technical Physics and Materials Science, 
H-1525 Budapest, P.O.Box ^9, Hungary 

The phase transitions of the recently introduced 2A 3A, 4A III reaction-difTusion model 
(G.Odor, PRE 69 036112 (2004)) are explored in two dimensions. This model exhibits site oc- 
cupation restriction and explicit diffusion of isolated particles. A reentrant phase diagram in the 
diffusion - creation rate space is confirmed in agreement with cluster mean-field and one-dimensional 
results. For strong diffusion a mean-field transition can be observed at zero branching rate charac- 
terized by a = 1/3 density decay exponent. In contrast with this for weak diffusion the effective 
2A 3A 4A reaction becomes relevant and the mean-field transition of the 2A 3A, 
2 A model characterized by q = 1/2 also appears for non-zero branching rates. 
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I. INTRODUCTION 

The classification of universality classes of nonequilib- 
rium systems is one of the most important tasks of sta- 
tistical physics [1,2]. Many of the known systems can 
be mapped onto some reaction-diffusion type of models, 
the behavior of them are the studied intensively in the 
past decades [3,4]. In these systems particle (A) cre- 
ation, annihilation and diffusion processes compete and 
by tuning the control parameters phase transition may 
occur from an active steady state to an inactive, absorb- 
ing state of zero density. For a long time only the critical 
"directed percolation" (DP) type of universality class has 
been known [5]. Later other classes were discovered re- 
lated to certain conservation laws or symmetries [6-9] , to 
long-range interactions [10-12], to boundary conditions 
[13-16] or to disorder [17-21]. These findings are all in 
agreement with the concepts of universality in equilib- 
rium systems. 

An extraordinary family of models has triggered a long 
debate among specialists recently [22-42]. The common 
behavior of these models is that for a particle production 
and annihilation at least two particles are needed (hence- 
forth they are called binary systems) and these reactions 
compete with the diffusion of isolated particles. Since for 
reactions at least a pair is needed while isolated particles 
can diffuse only these models can also be regarded as 
coupled systems [27]. The representative of this class is 
the so called diffusive pair contact process (PCPD) with 
reactions 2 A — > 3^1, 2 A [24]. The binary nature was 
found to be relevant in case of reactions of multi-species 
[33] too. 

The critical behavior of such models has been found to 
be different from all previously known classes (however 
there is still an ongoing debate on the precise values of 
critical exponents). The lack of symmetries, conserva- 
tion laws etc. have been motivating skepticism against 
the existence of a non-DP class transition and recently 



some studies suggested DP class behavior with extremely 
strong correction to scalings [44-46]. Field theoretical 
analysis [23] on the other hand indicate that the absence 
of the mass term corresponding to direct channel to the 
absorbing state {A 0) should be responsible for this 
"anomalous" behavior with respect to expectations based 
on equilibrium statistical physics. There is an other im- 
portant difference between binary systems and DP: there 
is no rapidity symmetry 



(1) 



between the field (0) and the response field (?/;) variables 
in the corresponding field theoretical description contrary 
the case of the DP process. Furthermore the lack of this 
relation is not the consequence of a symmetry breaking 
field of some boundary (like the i = boundary with 
long-ranged correlated order parameter field in case of 
pair contact process (PCP) [16]) or some disorder, but 
it is not there in the definition of these homogeneous, 
binary systems ^. 

Another odd feature is that bosonic (site unrestricted) 
and site restricted versions of these models show com- 
pletely different behavior. While site restricted mod- 
els investigated numerically exhibit the above continu- 
ous phase transition, the bosonic versions do not have 
steady state, but above an abrupt transition the density 
of particles diverges quickly [23,43]. The field theoreti- 
cal renormalization group (RG) analysis [23] predicts an 
upper critical dimension dc — 2, with logarithmic cor- 
rections at d = 2 for this class (PCPD). Simulations [32] 
have confirmed the mean- field scaling in two dimension in 
case of the 2A — 4^1, 2A — > binary production model. 

The site mean-field solution of general, 



nA A (n + k)A, mA {m ~ I) A, 



(2) 



models (with rt > 1, to > 1, A; > 0, ^ > and m — I > 0) 
resulted in a series of different universality classes de- 
pending on n and to [47]. This shows that above dc 



^Noh et al. claim that in their generalized PCPD model a long-range memory is generated by the diffusing isolated particles 
[41] 
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n and m are relevant parameters determining the 
type of continuous phase transitions. In particular for 
the n = m symmetrical case the density of particles above 
the critical point (ctc > 0) scales as 

p(x\(j-(Jcf, (3) 

with = 1, while at the critical point it decays as 

p oc i-" , (4) 

with a*^^ = = 1/n [38,47] (here "MF" de- 

notes mean-field value) . On the other hand for the n < m 
asymmetric case continuous phase transitions at zero 
branching rate ac = occur with 



In this work I continue the study of this model in 2d 
and show, that similar phase transition structure and 
critical behavior can be obtained. This is somewhat sur- 
prising, since one may expect that the diffusion is less 
relevant in higher dimensions due to its short interaction 
range. A very recent study using exact methods [49] 
showed that the particle density fluctuation and den- 
sity correlation function are diffusion dependent in the 
bosonic PCPD model for d > 2. In this work I give 
numerical evidence for diffusion dependence in a site re- 
stricted, binary model in d = 2. 



II. THE 2A 3A, 4A ^ MODEL 
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MF 
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(5) 



For n > m the mean-field solution provides first order 
transition. 

By going beyond site mean-field approximations it 
turns out that the above classification is not completely 
satisfying. In a previous paper [40] I investigated the 
2A — > 3A, AA — + model by cluster mean-field approxi- 
mations and simulations in Id and showed that the dif- 
fusion plays an important role: it introduces a different 
critical point besides the one at a = branching rate 
with eq. (5) exponents. The non-trivial critical point, 
obtained for low diffusion rate exhibits the universal be- 
havior of the 2A 3A, 2A ^ (PCPD) model owing 
to the generation of the effective 2A reaction via: 
2A ^ 3A 4:A ^ d) [48]. 




FIG. 1. Phase diagram of the 2A 3A ^ AA ^ 9 model. 
Stars correspond to = 2, boxes to A'' = 3, bullets to = 4 
and triangles to A'' = 5 cluster mean-field approximations. Di- 
amonds denote Id, -I- signs 2d simulation data, where PCPD 
class transitions are found. The lines serve to guide the eye. 
At the a — line asymmetric, eq. (5) type mean-field transi- 
tion occurs. 



This binary production reaction-diffusion model is de- 
fined by the rules: 



2A 3A 

4^^0 
v40 ^ 0A 



with rate a 
with rate X = I — a 
with rate D . 



(6) 



Here D denotes the diffusion probability and a is the pro- 
duction probability of the particles. The site occupancy 
is restricted to or 1 particle. In [40] the cluster mean- 
field approximations were determined on Id lattices for 
N = 1,2, ..5 cluster sizes. The corresponding reentrant 
phase diagram is shown on Fig. 1. Although cluster 
mean-field approximations based on d > 1 lattices may 
result in transition points at other locations, the univer- 
sal features are expected to be the same. Therefore I 
compare the simulation results with this approximation. 




FIG. 2. Simulation results for the steady state density at 
diffusions D = 0.5, 0.35, 0.1, 0.05 (solid lines from top to bot- 
tom) and N = 5 level cluster mean-field approximation data 
for D — 0.5, 0.05 (dashed lines from top to bottom). The 
insert shows the region near ct = magnified. 
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A. Simulation results 

I performed simulations in two dimensions in L = 
1 — 7 X 10'^ linear sized systems with periodic bound- 
ary conditions. The simulations were started from fully 
occupied lattices. One elementary Monte Carlo step con- 
sists of the following processes. A particle and a number 
xi e (0, 1) are selected randomly; if xi < D a site ex- 
change is attempted with one of the randomly selected 
empty nearest neighbors (nn). The time is updated by 
1/n, where n is the total number of particles. A parti- 
cle and a number X2 G (0,1) are selected randomly. If 
X2 < cr and if the number of nn particles is 1 or 2 or 
3, one new particle is created at an empty site selected 
randomly. If X2 > cr and the number of nn particles is 
greater than 2 four randomly selected neighboring parti- 
cles are removed. The time (t) is updated by 1/n again. 
The density of particles was followed up to t^ax ^ 10^ 
Monte Carlo steps (throughout the whole paper the time 
is measured by Monte Carlo steps (MCS)). 

As one can see on Fig. 2 simulation data and the 5- 
point approximations fit qualitatively well. In both cases 
for weak diffusion rates (for D <^ 0.1 in 2d simula- 
tions) reentrant phase transitions occur with CTc > 0, 
while for strong diffusions a single phase transition at 
(Tc = branching rate can be found. The transition lines 
of the cluster mean-field approximations do not converge 
towards the simulation line as in Id (see Fig.l), but the 
2d MC curve occurs at lower diffusions. But this is not 
surprising, since the cluster mean-field calculations are 
performed on Id lattices. 




FIG. 3. Density decay times t in the two-dimensional 
2A 3A, ^ model at D = 0.05. Different curves cor- 
respond to cr = 0.2715, 0.2708, 0.2704, 0.27, 0.2695, 0.269, 
0.2685, 0.268, 0.2677, 0.2675, 0.2673, 0.26715, 0.267, 0.2665, 
0.26 (top to bottom). The insert shows the corresponding 
local slopes. 



0.05 diffusion near the rightmost transition of Fig. 2 (at 
a ^ 0.27). By approaching Uc from the active phase the 
p{t)t^/'^ curves bend down rapidly for long times (beyond 
^ 10^ MCS). However this proved to be a finite size effect 
: the break-down of the density curves can be eliminated 
by increasing L. The largest system I could simulate had 
a linear size L = 7000. In this case no rapid and prema- 
ture curvatures was observed for t < 2 x 10^ MCS. As one 
can see on Fig. 3 for a > 0.2673 all curves veer up, while 
for cr < 0.2673 they veer down. A clear straight line - 
indicating scaling with the expected logarithmic correc- 
tion - can not be seen clearly. Even the a — 0.2673 curve 
shows some up and down curvatures in the last decade 
of the simulations. However as can be seen on the local 
slopes figure (see insert of Fig. 3) defined as: 



\n[p{t)/p{t/m)] 
\n(jn) 



(7) 



(where I used m = 2) the transition is around the ex- 
pected mean-field value of the PCPD class: a = 0.5 
[23,47]. Other curves exhibit strong curvatures for long 
times, i.e. for a > 0.2673 they veer up (active phase), 
while for a < 0.2673 they veer down (absorbing phase). 

The steady state density in the active phase near the 
critical phase transition point is expected to scale as 
p{oo) oc \(T — (Tc\'^ . Using the local slopes method one 
can get a precise estimate for f3 and see the corrections 
to scaling 



lnp(oo, aj) - lnp(oo, a^-i) 
In(CTi) - ln((Tj_i) 



(8) 



The steady state behavior at the Cc > transition for 
D — 0.05 was investigated using CTc = 0.2673(2) from 
the density decay analysis. Here the local slopes tend to 
(3eff — 0.98(2) without showing any relevant correction 
to scaling (see Fig. 4). This agrees with the mean-field 
value of the PCPD model again [23,47]. 
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FIG. 4. Peff as the function of a — CTc in the two dimen- 
sional 2A 3A, 4A model near the CTc = 0.2673 critical 
point for D — 0.05. The solid line shows a linear fitting. 



I explored the scaling behavior in more detail at D = 
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One may expect the same kind of transition all along 
the CTc > transition line. Indeed simulations showed 
that the density decays in a similar way at transitions 
with D = 0.01, 0.05, 0.09. 

To see the transition near ac — (horizontal axis on 
Fig.l) I determined the steady state value of p(oo, cr) for 
several as at D = 0.05 diffusion. The steady state den- 
sity was determined by running the simulations in the 
active phase near a = 0, by averaging over '--^ 100 sam- 
ples in a time window following the level-off is achieved. 
The smallest value I tested was a — 10^^, when I had 
to go up to t = 10^ MCS to reach a steady state (on 
a. L — 2000 sized system). By looking at the data it 
is quite obvious that the transition is at Cc = as the 
cluster mean-field approximations predicted. 

The effective order-parameter exponent (Fig. 5) tends 
to f3 = 0.505(5) as tJ ^ corroborating the cluster mean- 
field prediction: eq.(5). Assuming a correction to scaling 
of the form 

(3,ff^p-at-f^^ (9) 
fitting results in (3i — 0.5 as can be read-off from Fig. 5 . 



mean-field characteristics of the PCPD model because 
the effective 2A d) reaction (via 2A ^ 3A ^ 4A ^ 0) 
becomes relevant. The understanding of this diffusion 
dependence is a challenge for field theory. Similar reen- 
trant phase diagram has been observed in case of the 
unary production, triplet annihilation model (A 2A, 
3A 0) [51] and in a variant of the NEKIM model [52]. 
In all cases the diffusion competes with particle reaction 
processes, and the bare parameters should somehow form 
renormalized reaction rates which govern the evolution 
over long times and distances. An interesting question is 
whether this scenario extends above d = 2 dimensions as 
the cluster mean- field approximation predicts. A very re- 
cent non-perturbative RG study [53] finds similar phase 
diagram in case of the A 2A, 2A model for d > 3 
dimensions. That work points out that non-perturbative 
effects arise and there is a threshold {X/ D)th{d) above 
which DP, while below it a type (5) mean-field transition 
at = appears. 

The simulations also showed that at the CTc > tran- 
sition the finite size effects and corrections to scaling are 
very strong. I had to go up to 7000 x 7000 sized sys- 
tems and tjnax = 2 X 10^ MCS to see the appearance of 
the expected mean-field scaling with exponents a = 0.5, 
(3 = 1. Showing clear scaling for more than a decade with 
the predicted logarithmic corrections [23] is beyond the 
scope of this study, yet these simulation results for a 2d 
binary system are by far the largest scale ones published 
so far. On contrary the scaling at the Cc = critical 
point is clear with /3 = 0.505(5) and correction to scaling 
exponent /3i = 0.5. 
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FIG. 5. Peff as the function of a in the two dimensional 
2 A 3 A, 4 A model near the CTc = phase transition at 
D = 0.05. The solid line shows a linear fitting. 
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III. CONCLUSIONS 

In conclusion I have investigated the (D-a) phase di- 
agram of the two dimensional 2A 3A, 4A model 
with site restriction and explicit particle diffusion. Ex- 
tensive simulations gave numerical evidence that a reen- 
trant phase diagram emerges as in one dimension and 
predicted by cluster mean-field approximations [40] . This 
somewhat surprising results mean that diffusion plays rel- 
evant role even in d = 2 dimension. For high diffusion 
rates only a mean-field transition at cr = branching rate 
can be found, while for low diffusion an other transition 
type at CTc > appears. This latter transition shows the 



[1] U. Tauber, Adv. in Solid State Phys. 43 659 (2003). 

[2] G. Odor, Universality classes in nonequilibrium lattice 
systems, Rev. Mod. Phys. 76 (2004) in press, eprint: 
cond-mat/0205644. 

[3] For references see : J. Marro and R. Dickman, Nonequi- 
librium phase transitions in lattice models, Cambridge 
University Press, Cambridge, 1999. 

[4] H. Hinrichsen, Adv. Phys. 49, 815 (2000). 

[5] H. K. Janssen, Z. Phys. B 42, 151 (1981); P. Grassberger, 
Z. Phys. B 47, 365 (1982); G. Grinstein, Z-W. Lai and 
D. A. Browne, Phys. Rev A 40, 4820 (1989). 

[6] P. Grassberger, F. Krause and T. von der Twer, J. Phys. 
A:Math.Gen., 17, L105 (1984). 

[7] N. Menyhard, J.Phys.A 27, 6139 (1994); N. Menyhard 
and G. Odor, J. Phys. A. 28, 4505 (1995). 



4 



[8] J. L. Cardy and U. C. Tauber, J. Stat. Phys. 90, 1 (1998). 
[9] M. Rossi, R. Pastor-Satorras and A. Vespagnani, Phys. 
Rev. Lett. 85, 1803 (2000). 
[10] M. C. Marques and A. L. Ferreira, J. Phys. A 27, 3389 
(1994). 

[11] H. K. Janssen, K. Ocrding, F. van Wijland and H. Hil- 

horst, Eur. Phys. J. B 7, 137 (1999). 
[12] H. Hinrichsen and M. Howard, Eur. Phys. J. B 7, 635 

(1999). 

[13] P. Frojdh, M. Howard, K. B. Lauritsen, Int. J. Mod. 

Phys. B 15, 1761 (2001). 
[14] C. Kaiser and L. Turban, J. Phys. A 27, L579 (1994); 

C. Kaiser and L. Turban, J. Phys. A 28, 351 (1995). 
[15] H. Hinrichsen and G. Odor, Phys. Rev. E 58, 311 (1998). 
[16] I. Jensen and R. Dickman, Phys. Rev. E 48, 1710 (1993); 

I. Jensen, Phys. Rev. Lett. 70, 1465 (1993). 
[17] A. G. Moreira and R. Dickman, Phys. Rev. E 54, R3090 

(1996) . 

[18] A. J. Noest, Phys. Rev. Lett. 57, 90 (1986); A. J. Noest, 

Phys. Rev. B 38, 2715 (1988). 
[19] I. Wobman, D. ben Avraham, A. Cohen, and S. Havlin, 

Phil. Mag. B 77, 1401 (1998). 
[20] R. Cafiero, A. Gabrielli, and M. A. Munoz, Phys. Rev. E 

57, 5060 (1998). 
[21] J. Hooyberghs, F. Igloi and C. Vandorzande, Phys. Rev. 

Lett 90, 100601 (2003). 
[22] P. Grassberger, Z. Phys. B 47, 365 (1982). 
[23] M. J. Howard and U. C. Tauber, J. Phys. A 30, 7721 

(1997) . 

[24] E. Carlon, M. Hcnkol and U. Schollwock, Phys. Rev. E 

63, 036101-1 (2001). 
[25] H. Hinrichsen, Phys. Rev. E 63, 036102-1 (2001). 
[26] G. Odor , Phys. Rev. E62, R3027 (2000). 
[27] H. Hinrichsen, Physica A 291, 275-286 (2001). 
[28] G. Odor , Phys. Rev. E 63, 067104 (2001). 
[29] M. Henkel and U. Schollwock, J. Phys. A 34, 3333 (2001). 



[30] M. Hcnkcl and H. Hinrichsen, J. Phys. A 34, 1561 (2001). 
[31] K. Park, H. Hinrichsen, and In-mook Kim, Phys. Rev. E 

63, 065103(R) (2001). 
[32] G. Odor, M. A. Santos, M. C. Marques, Phys. Rev. E 

65, 056113 (2002). 
[33] G. Odor, Phys. Rev. E 65, 026121 (2002). 
[34] K. Park and I. Kim, Phys. Rev E 66, 027106 (2002). 
[35] G. Odor, Phys. Rev E 67, 016111 (2003). 
[36] J. Kockelkoren and H. Ghate, Phys. Rev. Lett. 90, 

125701 (2003). 

[37] R. Dickman and M. A. F. de Menenzes, Phys. Rev. E 66, 
045101 (2002). 

[38] K. Park, H. Hinrichsen and L Kim, Phys. Rev. E 66, 

025101 (2002). 
[39] G. Odor, Braz. J. of Phys. 33, 431 (2003). 
[40] G. Odor, cond-mat/0307351, Phys. Rev. E 69, 036112 

(2004). 

[41] J. D. Noh and H. Park, Phys. Rev. E 69, 016122 (2004). 

[42] For a recent review see: The non- equilibrium phase 
transition of the pair-contact process with diffusion, M. 
Henkel and H. Hinrichsen, cond-mat/0402433. 

[43] G. Odor and N. Menyhard, Physica D 168, 305 (2002). 

[44] H. Hinrichsen, Physica A 320, 249 (2003). 

[45] H. Hinrichsen, cond-mat/0302381. 

[46] G. T. Barkema and E. Carlon, cond-mat/0302151. 

[47] G. Odor, Phys. Rev. E 67, 056114 (2003). 

[48] This was confirmed by H. Ghate in a private communi- 
cation. 

[49] M. Paessens, G. M. Schuetz, cond-mat/0311568. 

[50] B. P. Lee, J. Phys. A 27, 2633 (1994). 

[51] R. Dickman, Phys. Rev. B 40, 7005 (1989); R. Dickman, 

Phys. Rev. A 42, 6985 (1990). 
[52] N. Menyhard and G. Odor, Phys. Rev. E 68, (2003) 

056106. 

[53] L. Canet, H. Ghate and B. Delamotte, cond- 
mat/0403423. 



5 



